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Abstract 


Work is currently underway at the NASA Lewis Research Center to develop an analytical method for 
predicting the performance degradation of a helicopter operating in icing conditions. A brief survey is performed 
of possibilities available to perform such a calculation along with the reasons for choosing the present approach. 
A complete description of the proposed jobstream is given as well as a discussion of the present state of the 
development. 
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Damping Length Constant 
Interaction Coefficient Matrix 
Local Ice Thickness, m 
Total Collection Efficiency 
Freezing Fraction 

Dimensionless Sand Grain Roughness 
Equivalent Sand Grain Roughness, m 
Mixing Length 

Mass Flow Rate Due To Impinging Liquid, kg/sec 

Mass Flow Rate Which Freezes, kg/sec 

Mass Flow Rate Due To Runback, kg/sec 

Rotor Blade Radius, m 

Static Temperature, °C 

Inviscid Velocity, m/sec 

Tangential Edge Velocity, m/sec 

Total Normal Velocity Component, m/sec 

Total Tangential Velocity Component, m/sec 

Friction Velocity, m/sec 

Normal Velocity, m/sec (Boundary Layer Calculations) 

Freest ream Velocity, m/sec 

Wake Induced Normal Velocity, m/sec 

Wake Induced Tangential Velocity, m/sec 

Distance Along Surface, m 

Distance Normal To Surface, m 

A Reynolds Number, y i^/V 

Total Flow Angle, deg. 

Shaft Angle, deg. 

Local Collection Efficiency 
Flapping Angle, deg. 

Airfoil Segment Length, m 
Icing Time, sec 
Ice Density, g/m 3 


6* Displacement Thickness, m 

5u c Perturbation Velocity Due To Viscous Effects, m/sec 

0 Kinematic Angle of Attack, deg. 

k Universal Constant, Also Sweep Parameter 

u Kinematic Viscosity, m 2 /sec 

<f> Aerodynamic Flow Angle, deg. 

^ Azimuth Angle, deg. 

0 Angular Velocity, rad/sec 


Introduction 

Historically, certification/qualification of a helicopter for flight into known icing conditions has been a 
problem. This is because of the current emphasis on flight testing for verification of system performance. Flight 
testing in icing conditions is difficult because, in addition to being dangerous and expensive, many times conditions 
which are sought after cannot be readily found in nature. The problem is compounded for helicopters because of 
their small range in comparison to many fixed wing aircraft. Thus, helicopters are forced to wait for conditions 
to occur in a certain region rather than seeking them out. These and other drawbacks to flight testing have 
prompted extreme interest in developing validated alternatives to flight testing. One such alternative is theoretical 
prediction. 

Any complete analysis of this problem must be able to perform three calculations: 

1) Clean performance of helicopter main rotor; 

2) Characterization of the effect of icing on performance; 

3) Prediction of ice shedding. 

The ability to predict the clean performance of a helicopter main rotor has existed for some time. Traditionally 
lifting line analyses have been the main vehicle for this calculation but recently Navier-Stokes analyses have become 
a potential option. 

The areas of difficulty lie in the last two calculations: prediction of the effect of icing on performance and 
prediction of ice shedding. Prediction of the effect of icing on rotor performance is possible using empirically based 
correlation methods. Until recently, this has been the primary means of determining the effect of icing on rotor 
lift, drag, and moment characteristics. The empirical methods have the advantage of being very simple and yielding 
acceptable results for their range of applicability. The main drawback of these methods however, is that they are 
generally limited in terms of the conditions under which they can be properly applied. Thus, a valid analytical 
method would have a distinct advantage over the correlations. Two possible analytical tools for determining the 
lift, drag, and moment change due to ice accretion are the Interactive Boundary Layer 1 (IBL) procedure and Navier 
Stokes 2 analysis. As with any analytical method both of these require characterization of ice shape in order to 
perform calculations on the effective airfoil shape. While showing a great deal of promise, the Navier Stokes 
analyses are all still computationally intensive. The IBL procedure seems appropriate for the task because it is not 
computationally intensive and requires no grid. Thus, it was decided to use the IBL procedure for the current work. 

A great deal of research into the ice shedding phenomena is currently underway. It is hoped that eventually 
the capability will exist to analytically predict the time and location of shedding events on a rotor. Currently 
however, prediction schemes still rely heavily on correlation methods. 

Proposed Method of Analysis 

As discussed earlier, prediction of the effect of icing on the performance of a helicopter main rotor is a 
complicated task requiring several steps of calculation. Figure 1 shows a flowchart which outlines the various steps. 
First, the clean performance of the main rotor must be determined as a baseline. Because other codes which will 
be used in this analysis are two dimensional, a performance code which makes use of blade element theory allows 
for easy transfer of information. This type of performance code relies on experimental data tables for two 
dimensional airfoil lift, drag, and moment values which are integrated to obtain the overall rotor performance. Once 
the clean performance has been determined, the trimmed values of Mach number and angle of attack can be fed into 
LEWICE 3 , the accretion analysis. LEWICE uses an azimuthal average of these values and calculates the ice shapes 


2 


ftl various radial locations along the blade. This averaging technique is discussed in more detail in a later section. 
The new "effective" airfoil shapes are then passed to the IBL procedure which calculates the new lift, drag, and 
moment characteristics. A check is then made to determine if any shedding has occurred and if so, where. If 
shedding has occurred, it is assumed that the radial location of the event is now "clean" and has normal lift, drag, 
and moment values. It should be noted that the process of ice shedding often leaves residual ice still attached to 
the surface. Currently, no method exists for predicting how much residual ice will remain after a shedding event. 
The new characteristics are then used by the performance code instead of the original data tables to calculate the 
"degraded" performance. Although not included in the present study, future efforts will attempt to predict the 
effects of various deicing/anti-icing systems using codes now under development. 4 The following sections discuss 
in more detail the calculation procedure for each step. 

Helicopter Performance Prediction 

The central part of the proposed jobstream is the performance prediction code. Although other methods 
exist, it has been decided to use Boeing Helicopters’ B65 performance code, which uses lifting line theory in its 
aerodynamic model. Because the code is proprietary in nature, discussion will be limited to a brief description of 
lifting line theory in general. 

Lifting line theory makes use of blade element theory in which the rotor blade is broken into several 
discrete radial sections. Characteristics are calculated for each section as a function of radial and azimuth location 
and then integrated to obtain rotor performance values. Each section acts as a steady 2-D airfoil section with three 
dimensional, compressibility, and unsteady effects, including dynamic stall, usually included as corrections to the 
overall behavior. Section characteristics are normally obtained as a function of Mach number and angle of attack 
by use of an empirical set of data tables. 

The crux of a lifting line performance code is the calculation of the flow angle of attack over a 2-D slice 
of the rotor blade. This flow angle is given by 


where <j> arises from the purely aerodynamic behavior and 0 comes from the blade aeroelasticity. The aeroelastic 
behavior of the rotor is usually significant because the high aspect ratio blades tend to be very flexible and respond 
to the various airloads. A typical velocity triangle for a rotor section is given in Figure 2. Here the total normal 
velocity component, u^, and the tangential component, %, are given by 

d$ f 

u p =V cos(a>inf (2) 

« r = QR + Ksinta^^ilr + V T 


In a normal calculation the rotational speed, OR, and shaft angle, a g , are known. Thus, the remaining quantities 
which need to be calculated are the flapping velocity and the normal and tangential components of the wake induced 
velocity. The values are obtained through various procedures and iterated upon until the desired trimmed thrust 
level is achieved. A more complete description of blade element theory can be found in Reference 5 and an in depth 
discussion of a typical lifting line analysis can be found in Reference 6. 

Ice Accretion Prediction 

A numerical analysis LEWICE 3 has been developed by the NASA Lewis Research Center which has the 
ability to predict the analytical ice shape which accretes on a given component exposed to an icing condition for a 
known period of time. This analysis models four critical steps in the icing process, which are: 

(1) Flowfield calculation about the component; 

(2) Water droplet impingement characteristics; 

(3) Heat transfer processes; 
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(4) Ice accumulation normal to the surface. 

Each of these areas are discussed briefly in the following sections. 

Flowfield Calculation 

Predicting the flowfield about a body which has an accreted ice shape presents several challenges because 
of the irregular effective airfoil shape which often occurs. The potential flow program developed by Hess and 
Smith 7 is incorporated into LEWICE. This method makes use of distributed sources, sinks, and/or vortices to 
describe the flowfield about a body which has been modeled by a series of line segments. Comparisons of results 
to experimental clean body data have been favorable for the normal ranges of incompressible flow. A potential 
problem does exist for applying the assumption of incompressible flow to regions near the tip of a rotor blade. It 
is felt that the high centrifugal forces in this region will cause shedding to occur, especially for wanner 
temperatures. Thus, this particular limitation is not seen as causing any significant difficulty. 

Impingement Characteristics 

Of primary importance in any ice accretion analysis is characterization of the region of impinging water 
droplets. This characterization consists of the limits of impingement as well as the distribution of the mass of the 
impinging liquid. This is typically obtained by performing an analysis of droplet trajectories from far upstream. 
Two important parameters which arise are the total and local collection efficiency. The total collection efficiency 
is defined as the ratio of the total mass of impinging liquid over the theoretical mass of impinging liquid which 
would occur if all of the droplet trajectories were straight lines. The local collection efficiency is based on the same 
definition, except that it pertains to a specific location on the body. A pictorial representation of this is given in 
Figure 3. 

The droplet trajectory analysis used in LEWICE is based on the work of Frost, Chang, Shieh, and Kimble. 8 
The method has the ability to calculate trajectories and impingement characteristics of an arbitrarily shaped particle. 
Although generality has been maintained, for this application it can usually be assumed that the particles are 
spherical and gravity forces are negligible. This simplifies the calculation somewhat. 

Heat Transfer Processes 

The freezing process is modeled in LEWICE by performing a mass and energy balance on a control volume 
located on the surface and extending beyond the boundary layer, as shown in Figure 4. Each segment which 
describes the surface of the accreting component has a corresponding control volume. The runback model first 
developed by Messinger 9 is incorporated here. An important quantity in this analysis is the freezing fraction which 
is the ratio of freezing liquid within a control volume over the total amount of liquid entering. The freezing fraction 
can be expressed as: 


/= 


m, 


m+m_ 


( 3 ) 


Once the temperature and the freezing fraction are known, the mass balance calculates the mass flow of the liquid 
runback out of the control volume. Any liquid which leaves the control volume is assumed to leave in the direction 
away from the stagnation point. Surface roughness strongly influences the local heat transfer processes. The 
equivalent roughness concept is used here. This concept models the actual surface roughness by using an average 
value which yields the same heat transfer characteristics. This aspect of the analysis is considered a weak point and 
much effort is being expended to improve the current heat transfer model in general. 10 

Surface Ice Growth 

The ice is assumed to grow normal to the surface. The ice growth rate can be defined using the expression 
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for the freezing fraction and is given as: 




(4) 


This can be redefined to yield an expression for the ice thickness as: 
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Pi 


( 5 ) 


The iced component is then calculated by adding the corresponding ice thickness normal to each matching segment. 
This results in a new airfoil surface coinciding with the specified icing time. 

Calculation Procedure 

Using LEWICE to predict a 2-D ice accretion at a specified radial location on a helicopter rotor is not a 
straightforward calculation. Thus, a brief discussion of the procedure used in this calculation is warranted. 
LEWICE is a steady state code designed for use on a fixed wing where the velocity and angle of attack are constant 
throughout a given time step. However, for application to a helicopter in forward flight where the local angle of 
attack and Mach number are constantly changing some averaging procedure is necessary. A technique developed 
in 1983 by Korkan, Dadone, and Shaw 11 dealt with this problem. Here, while attempting to simplify the analysis 
of a helicopter main rotor in forward flight with a rime ice accretion, several methods of averaging were 
investigated. It was found that if the local Mach number and angle of attack produced by a helicopter performance 
code were averaged and input into an icing performance degradation analysis the predicted change only differed by 
±2% over that of the traditional method of calculating values at specific azimuth locations around the disk. In view 
of this, a similar averaging technique is employed in the present study. First, trimmed performance values are 
calculated using a helicopter performance code. The local flow angle at the desired radial location required for trim 
is then averaged azimuthally. The local velocity is taken to be the rotational velocity at the specified radial location. 
This is, in effect, the averaged velocity. These values are then used in LEWICE as the velocity and angle of attack. 
The calculation procedure from then on is carried out as any normal LEWICE calculation using established 
guidelines given in the LEWICE User’s Manual. 3 As shown in Figures 5 and 6, previous comparisons of results 
obtained using this procedure to experiment have been very good. More in depth comparisons using this procedure 
are given in Reference 12. 

Performance Penalty Prediction 


Critical to predicting the performance degradation of an icing encounter on an aircraft is the ability to 
accurately compute the associated changes in sectional lift, drag, and moment characteristics. Navier-Stokes 
methods have this ability. 2 However, current Navier-Stokes schemes require a significant amount of computer time. 
A great deal of research is being performed in this area and as mainframe technology advances and the efficiency 
of the schemes improve it is anticipated that Navier-Stokes methods will eventually be the means for this type of 
calculation. However, at the present, a simpler short term solution is needed. The Interactive Boundary Layer 
(IBL) procedure developed by Cebeci 1 has been chosen for its simplicity and apparent success at predicting drag 
values. The IBL procedure is attractive in that it does not require a computational grid. 

The IBL approach consists of solution of inviscid flow and boundary layer equations which are coupled 
so that one influences the other. The inviscid calculations are derived from the panel method developed by Hess 
and Smith. 7 This method models the airfoil and ice shape by a series of line segments which contain distributed 
source and/or vorticity strength. The set of simultaneous linear equations are solved such that the normal velocity 
boundary condition at the midpoint of the segments is satisfied. The total normal velocity at each segment midpoint 
is zero for inviscid flow. When modeling the effects of the boundary layer however, the normal velocity, v n , is 
given by the derivative along the surface of the product of the displacement thickness and the tangential velocity, 
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expressed as: 


d(u € V) 


ds 


( 6 ) 


This surface blowing distribution has the effect of displacing the dividing streamline outward a distance equal to the 
displacement thickness. Thus, this approach is equivalent to the classic procedure of maintaining the zero normal 
velocity boundary condition but redefining the aerodynamic surface to include the displacement thickness. While 
the Kutta condition still must be maintained, researchers have found that better results are obtained if it is applied 
to the displacement surface rather than the original surface. 

The boundary layer equations for steady two-dimensional incompressible flows are solved with the velocity 
distribution at the edge of the boundary layer coming from inviscid flow theory. The total velocity distribution can 
be expressed in terms of the inviscid velocity and the perturbation velocity due to viscous effects: 

«,(*>=« '(*)+&«.(*) w 


where: 



( 8 ) 


The range x # x <* x b is normally taken to be the airfoil plus two chord lengths downstream. Equation 8 
provides an outer boundary condition for the viscous flow calculations and represents the interaction between the 
viscous and inviscid flow. Equation 8 can be generalized to the form: 


«,(*)=«:<*)+£ cjL(u t y) r (u t i>-)]\ 


(9) 


where u e *(x) is the inviscid velocity distribution containing the displacement thickness effect calculated from the 
previous sweep. The interaction coefficient matrix, Cy, is obtained from a discrete approximation to the Hilbert 
integral. 

The boundary layer solution procedure is derived from Keller’s box scheme. The second order finite 
difference approximations are written in terms of Falkner-Skan variables. Solutions with separation are computed 
using the inverse form of the equations. The FLARE approximation of Reyhner and FJugge-Lotz 13 is used which 
sets the convective term, u(du/3x), in the recirculation region to zero. This eliminates the numerical instabilities 
associated with integrating the boundary layer equation against the direction of local flow. The inaccuracies 
resulting from this approximation are generally considered small because magnitude of the values of u in the 
reversed flow region are small compared to the external flow velocity. Although methods exist for improving the 
numerical method if necessary, no attempt was made to do so here. The finite difference approximations yield a 
nonlinear system of algebraic equations which are linearized by Newton’s method and solved by a block elimination 
procedure. This procedure is described in more detail in Reference 14. 

In order to deal with surface roughness associated with ice the mixing length expression of the Cebeci- 
Smith model 13 has been modified as 


L=K(y+Ay){l 



(10) 
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where Ay is a function of the equivalent sand-grain roughness k,. Ay can be expressed in terms of dimensionless 
quantities. 


Ay *=0.9 
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and 



( 12 ) 


(13) 


The roughness is converted into equivalent sand-grain roughness by using the procedure of Smith and Kaups 16 and 
assuming the ratio of equivalent sand-grain roughness to the roughness of the applied elements to be a function of 
the concentration and shape of the roughness elements. More detail on this procedure is given in Reference 1 . 

Accreted ice on an airfoil can substantially alter the leading edge geometry in a short period of time which 
causes rapid variations in flow properties. Thus, some difficulty is met trying to obtain acceptable solutions from 
the inviscid and viscous flow calculations. Thus, "blanketing" and continuation techniques have been employed to 
minimize this problem. A detailed description of these techniques is given in References 1 and 16. 

Work by Shin, et at 1 has shown that the IBL procedure can acceptably predict drag values of iced airfoils. 
Results from that work are shown in Figures 7, 8, and 9. It was found that the IBL procedure had a tendency to 
underpredict the experimental drag values. This was due in large part to the fact that the theoretical ice shape from 
which the IBL procedure made its calculations was much smoother than the experimental shape. In any scheme 
of this nature the quality of the drag predictions will necessarily depend upon the ice shape prediction and the ability 
to account for roughness. 

Ice Shedding Prediction 

An icing analysis of a rotating system differs from that of a fixed system in that ice shedding becomes a 
predominant factor. The combination of centrifugal force and vibratory airloads makes shedding commonplace for 
a helicopter main rotor. In a general sense, ice shedding occurs when the centrifugal, bending, vibratory, and 
aerodynamic forces acting on a mass of ice causes the stress within the ice to exceed a critical value. When this 
critical value is surpassed, failure occurs within the ice and aerodynamic forces carry it away. Since LEWICE 
provides the ability to predict the ice shape and approximations exist for calculating the ice density, estimation of 
the centrifugal force acting on a section of ice is possible. Calculation of the stresses resulting from this force can 
be carried out with levels of complexity ranging from a simplistic approximation to a complete Finite Element 
Analysis. Considerations such as CPU time and FORTRAN compatibility have led to the decision to use an 
approximate model. In the present approach, only centrifugal forces are considered in the stress calculations. The 
stresses resulting from bending, vibratory, and aerodynamic loads are assumed to be small. Ignoring the bending 
loads is appropriate for this phase of the study because the initial comparisons will be made to powered force model 
data taken in the NASA Lewis Research Center Icing Research Tunnel. The blades used for this experiment were 
extremely stiff and it is anticipated that very little bending of the blades actually occurred. Recent work by 
Scavuzzo 18 has shown that aerodynamic loading can be significant enough for consideration. Future efforts will 
attempt to include the aerodynamic forces in the overall calculation of shear and normal stresses. Here, the normal 
and shear stresses are computed based on the ice shape areas provided by LEWICE. Calculations begin at the tip 
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and progress inward towards the hub. Once the stress levels have exceeded the failure stress, it is assumed that 
all of the ice from that location to the tip sheds away. The calculation scheme allows ice to begin accreting again 
in the region after the shed. 

Determination of the failure stress of the accreted ice is by far the more difficult task. Scavuzzo, et at 8 9 
has attempted to experimentally determine the critical shear and normal stresses of accreted ice for various 
conditions. These results show a strong dependence on surface temperature (above -11 °C) and surface roughness 
as well as some dependence on wind velocity and droplet size. Historically, as is the case here, these types of 
experiments have shown a great deal of scatter in the data. Thus, no method currently exists which will absolutely 
determine when and where a shedding event will occur. Currently, only probabilities of shedding events can be 
calculated. For the purposes of this study if the probability of shedding is greater than 50% (based on Scavuzzo’s 
data), then the ice is assumed to have shed. 


Summary 

A method of computing the effects of an icing encounter on the performance of a helicopter main rotor has 
been proposed. This method makes use of several codes, each of which play a vital part in the overall calculation. 
The Boeing Helicopters performance code, B65 is used to make the actual performance calculations. NASA Lewis’s 
ice accretion code, LEWICE is used to calculate the shape of the accreted ice at various radial locations along the 
rotor. The IBL procedure of Cebeci determines the effect of the ice accretion on the lift, drag, and moment 
characteristics of two dimensional airfoil slices along the radius. Shedding events are predicted by approximating 
the stresses in the ice along the blade and comparing to experimental data obtained by Scavuzzo. All of the codes, 
with the exception of the shedding model, have been independently compared with experiment with generally 
favorable results. Currently, the performance code has been linked with the ice accretion code. Efforts are now 
underway to complete the jobstream by linking the IBL procedure and the shedding model to the existing 
performance/accretion code. Once completed, the results will be compared to experimental data obtained in the 
NASA Lewis Icing Research Tunnel. 
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Figure 1.— Helicopter performance in icing prediction jobstream. 



Figure 2 —Blade section velocity triangle. 
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S u = Upper - surface impingement limit 
= Lower - surface impingement limit 
H = Forward projection of the airfoil height 

Total collection efficiency 



Total collection efficiency „ _ 

P = ds 

Figure 3. — Definition of total and local collection efficiency (Ref. 3). 
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Figure 5. — Comparison between experiment and theoretical 
prediction of LEWICE. (Ref. 9) 
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Figure 6. — Comparison between experiment and theoretical 
prediction of LEWICE. (Ref. 9) 
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Figure 7. — Comparison between experiment and theoretical 
prediction of 1BL procedure (Ref. 15). (Line is a fit for 
experimental data). 



Figure 8. — Comparison between experiment and theoretical 
prediction of 1BL procedure. V = 58 m/sec, T s = -27.8 °C, 
LWC = 1 .0 g/m3, MVD = 1 2^ m, t = 5 min) (Ref. 1 5). 
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Figure 9. — Comparison between experiment and theoretical 
prediction of IBL procedure. (V = 58 m/sec, T s = -9.67 
°C, LWC = 1 .3 g/m3, MVD = 20^ m,x = 5 min.) (Ref. 1 5) 
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